Test SEIR Model

Testing data

fips_table <- get_fips(datadir = datadir)
reopen.df <- get_reopendata(datadir=datadir)
alldat.df <- get_testdata(mindate = "2020-03-17")

reopen.df <- get_reopendata(datadir = datadir)
alldat.df$valueperM <- alldat.df$value*1e6/
  fips_table[alldat.df$state,"pop"]
allcumdat.df <- get_cumul_testdata(mindate = "2020-03-17")
allcumdat.df$valueperM <- allcumdat.df$value*1e6/
  fips_table[allcumdat.df$state,"pop"]
dat.df <- subset(alldat.df,state==statenow)
cumdat.df <- subset(allcumdat.df,state==statenow)

ggplot(dat.df)+geom_point(aes(x=numDate,y=value,color=variable))+
  scale_y_log10()+
  geom_vline(aes(linetype=ReopenType,xintercept=numDate),
             data=subset(reopen.df,State.Abbr==statenow))+
  ggtitle(paste(statenow,"Daily"))

ggplot(cumdat.df)+geom_point(aes(x=numDate,y=value,color=variable))+
  scale_y_log10()+
  geom_vline(aes(linetype=ReopenType,xintercept=numDate),
             data=subset(reopen.df,State.Abbr==statenow))+
  ggtitle(paste(statenow,"Cumulative"))

reopen.df$state <- reopen.df$State.Abbr
ggplot(alldat.df)+geom_point(aes(x=numDate,y=value,color=variable))+
  scale_y_log10()+facet_wrap(~state)+
  xlim(range(alldat.df$numDate))+
  geom_vline(aes(linetype=ReopenType,xintercept=numDate),
             data=reopen.df)+
  ggtitle("Daily")

ggplot(allcumdat.df)+geom_point(aes(x=numDate,y=value,color=variable))+
  scale_y_log10()+facet_wrap(~state)+
  xlim(range(alldat.df$numDate))+
  geom_vline(aes(linetype=ReopenType,xintercept=numDate),
             data=reopen.df)+
  ggtitle("Cumulative")

ggplot(alldat.df)+geom_point(aes(x=numDate,y=valueperM,color=variable))+
  xlim(range(alldat.df$numDate))+
  geom_vline(aes(linetype=ReopenType,xintercept=numDate),
             data=reopen.df)+
  scale_y_log10()+facet_wrap(~state)+
  ggtitle("Daily")

ggplot(allcumdat.df)+geom_point(aes(x=numDate,y=valueperM,color=variable))+
  xlim(range(alldat.df$numDate))+
  geom_vline(aes(linetype=ReopenType,xintercept=numDate),
             data=reopen.df)+
  scale_y_log10()+facet_wrap(~state)+
  ggtitle("Cumulative")

MCSim simulation

Set up MCSim file.

mdir <- "../MCSim"
source(file.path(mdir,"setup_MCSim.R"))
# Make mod.exe (used to create mcsim executable from model file)
makemod(mdir) 

## The mod.exe had been created.

model_file<- "SEIR.reopen.model.R"
exe_file<-makemcsim(model_file,modeldir=modeldir,mdir="../MCSim")

## * Created executable program '../model/mcsim.SEIR.reopen.model.R.exe'.

Deterministic Model

Start time = 60. Texas population.

set.seed(314159)
in_file <- "SEIR_Testing.in.R" 
out_dat <- mcsim(exe_file = exe_file,
                        in_file = in_file,
                 resultsdir = resultsdir)

## Execute: ./../model/mcsim.SEIR.reopen.model.R.exe ./SEIR_Testing.in.R ./sim.out

out_dat <- data.table(out_dat)
states_S.df <- melt(out_dat[,c("Time","S","S_C","Tot")],id.vars=1)
ggplot(states_S.df,aes(x=Time,y=value,color=variable))+geom_line() + scale_y_log10()+facet_wrap(~variable)+theme(legend.position = "none")

## Warning: Transformation introduced infinite values in continuous y-axis

states_E.df <- melt(out_dat[,c("Time","E","E_C")],id.vars=1)
ggplot(states_E.df,aes(x=Time,y=value,color=variable))+geom_line() + scale_y_log10()+facet_wrap(~variable)+theme(legend.position = "none")

## Warning: Transformation introduced infinite values in continuous y-axis

states_I.df <- melt(out_dat[,c("Time","I_U","I_C","I_T")],id.vars=1)
ggplot(states_I.df,aes(x=Time,y=value,color=variable))+geom_line() + scale_y_log10()+facet_wrap(~variable)+theme(legend.position = "none")

## Warning: Transformation introduced infinite values in continuous y-axis

states_R.df <- melt(out_dat[,c("Time","R_U","R_T","F_T")],id.vars=1)
ggplot(states_R.df,aes(x=Time,y=value,color=variable))+geom_line() + scale_y_log10()+facet_wrap(~variable)+theme(legend.position = "none")

## Warning: Transformation introduced infinite values in continuous y-axis

totals.df <- melt(out_dat[,c("Time","CumInfected", "CumPosTest", "CumDeath")],id.vars=1)
ggplot(totals.df,aes(x=Time,y=value,color=variable))+geom_line() + scale_y_log10()+facet_wrap(~variable)+theme(legend.position = "none")

## Warning: Transformation introduced infinite values in continuous y-axis

timedep.df <- melt(out_dat[,c("Time","ThetaFit","HygieneFit","FTraced", "lambda", "lambda_C", "rho_C", "delta", "c", "beta","Rt","Refft")],id.vars=1)
ggplot(timedep.df,aes(x=Time,y=value,color=variable))+geom_line() +facet_wrap(~variable,scales="free_y")+theme(legend.position = "none")

const.df <- melt(out_dat[,c("Time","NInit", "TIsolation", "R0", "c0", "TLatent", "TRecover", "IFR", "TStartTesting", "TauTesting", "TTestingRate", "TContactsTestingRate", "TestingCoverage", "TestSensitivity", "ThetaMin", "TauTheta", "PwrTheta", "HygienePwr", "FTraced0", "TPosTest", "TFatalDeath", "alpha", "kappa", "rho", "lambda0", "lambda0_C", "rho0_C","beta0")],id.vars=1)

## Warning in melt.data.table(out_dat[, c("Time", "NInit", "TIsolation", "R0", :
## 'measure.vars' [NInit, TIsolation, R0, c0, ...] are not all of the same type. By
## order of hierarchy, the molten data value column will be of type 'double'. All
## measure variables not of type 'double' will be coerced too. Check DETAILS in ?
## melt.data.table for more on coercion.

ggplot(const.df,aes(x=Time,y=value,color=variable))+geom_line()+scale_y_log10()+facet_wrap(~variable)+theme(legend.position = "none")

out_dat.tmp<-out_dat[,c("Time","N_pos","D_pos")]
names(out_dat.tmp)<-c("Time","positiveIncrease","deathIncrease")
obs.df <- melt(out_dat.tmp,id.vars=1)
ggplot(obs.df,aes(x=Time,y=value,color=variable))+geom_line() + 
  scale_y_log10(limits=c(1,NA)) +
  geom_point(data=dat.df,aes(x=numDate,y=value,color=variable))

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 23 rows containing missing values (geom_path).

## Monte Carlo test

Counts how many integration failures occur using 5000 random parameter sets.

# Without integration
in_file0 <- "SEIR_Testing_MTC0.in.R" 
out_dat0 <- mcsim(exe_file = exe_file,
                        in_file = in_file0,
                  out_file = file.path(resultsdir,"simMTC0.out"),
                 resultsdir = resultsdir,ignore.stdout = FALSE)
out_dat0 <- fread(file.path(resultsdir,"simMTC0.out"),select=1:25)
# out_mtc0.df <- melt(out_dat0,id.vars=1)
# out_mtc0.df$dointeg <- FALSE

# With integration
in_file <- "SEIR_Testing_MTC.in.R" 
out_dat <- mcsim(exe_file = exe_file,
                        in_file = in_file,
                  out_file = file.path(resultsdir,"simMTC.out"),
                 resultsdir = resultsdir,ignore.stdout = FALSE)
out_dat <- fread(file.path(resultsdir,"simMTC.out"),select=1:25)

# ggpairs(log10(out_dat[,2:17]),
#         lower = list(continuous = wrap("points", size=0.2,alpha=0.4)),
#         title="Succeeded integration"
#         )

integdat<-setdiff(out_dat0,out_dat)
# ggpairs(log10(integdat[,2:17]),
#         lower = list(continuous = wrap("points", size=0.2,alpha=0.4)),
#         title="Failed integration"
#         )
nfail <- nrow(integdat)
print(nfail)

## [1] 0

if (nfail > 50) {
  integdat$integ <- TRUE
  tmp <- out_dat
  tmp$integ <- FALSE
  all_mtc <- rbind(integdat,tmp)
  
  all_mtc[,2:22]<-log10(all_mtc[,2:22])
  ggpairs(all_mtc[,2:23],
          mapping = ggplot2::aes(color = integ,alpha=integ),
          lower = list(continuous = wrap("points", size=0.2))
          )
}
system(paste("rm",file.path(resultsdir,"simMTC.out")))

Test run using generic prior (used in validation runs)

resultsdir <- "TX.val"
fips_table <- get_fips(datadir = datadir)
dat.df <- get_testdata()
cumdat.df <- get_cumul_testdata()
burnin <- 0.1
datadatemax <- "2020-04-30"
priorfile<-"SEIR.reopen_priors_MCMC.in.R"
statenow <- "TX"
popnow <- fips_table$pop[fips_table$Alpha.code==statenow]
prior_template <- make_infile_template(dat.df,
                                      fips_table,
                                      state_abbr=statenow,
                                      usestatename = TRUE,
                                      createdir = TRUE,
                                      pathdir = resultsdir,
                                      priordir = "../priors",
                                      priorfile = gsub("MCMC","MTC",priorfile),
                                      usemobility = FALSE,
                                      mobilitydir = "../MobilityMetrics",
                                      isprior = TRUE
)
prior_dat <- mcsim(exe_file = exe_file,
                      in_file = prior_template,
                 out_file = file.path(resultsdir,gsub(".in.R",".out",prior_template)),
               resultsdir = resultsdir)

## Execute: ./../model/mcsim.SEIR.reopen.model.R.exe TX.val/SEIR_TX_MTC.in.R TX.val/SEIR_TX_MTC.out

  infile_template <- make_infile_template(dat.df,
                                        fips_table,
                                        state_abbr=statenow,
                                        usestatename = TRUE,
                                        createdir = TRUE,
                                        pathdir = resultsdir,
                                        X_iter = "2000",
                                        X_print = "10",
                                        priordir = "../priors",
                                        priorfile = priorfile,
                                        datadatemax=datadatemax,
                                        usemobility = FALSE,
                                        mobilitydir = "../MobilityMetrics")
  set.seed(exp(2))
  out_dat <- mcsim(exe_file = exe_file,
                          in_file = infile_template,
                   resultsdir = resultsdir)#,ignore.stdout = TRUE)

## Execute MCMC: ./../model/mcsim.SEIR.reopen.model.R.exe TX.val/SEIR_TX_MCMC.in.R TX.val/SEIR_TX_MCMC1.out

## * Created 'TX.val/SEIR_TX_MCMC1.out'

## Execute MCMC Check: ./../model/mcsim.SEIR.reopen.model.R.exe TX.val/SEIR_TX_MCMC.in.R TX.val/SEIR_TX_MCMC1.out

## * Created 'TX.val/SEIR_TX_MCMC1.check.out' from the last iteration.

  make_diagnostic(out_dat, subset(dat.df,state==statenow), 
                  subset(cumdat.df,state==statenow), burnin=0.1,
                  pdfname=file.path(resultsdir,
                                    paste0("Test.Validation.",
                                           statenow,".pdf")))

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 6 rows containing missing values (geom_col).

## Warning: Removed 30 rows containing missing values (geom_path).

## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 1935 rows containing missing values (geom_point).

## quartz_off_screen 
##                 2

Test run using state-specific priors (used for prediction runs)

statenow <- "TX"
resultsdir <- "TX.pred"
fips_table <- get_fips(datadir = datadir)
dat.df <- get_testdata()
cumdat.df <- get_cumul_testdata()
burnin <- 0.1
popnow <- fips_table$pop[fips_table$Alpha.code==statenow]
priorfile<-"SEIR.reopen_state_priors_MCMC.in.R"
prior_template <- make_infile_template(dat.df,
                                        fips_table,
                                        state_abbr=statenow,
                                        usestatename = TRUE,
                                        createdir = TRUE,
                                        pathdir = resultsdir,
                                        priordir = "../priors",
                                        priorfile = gsub("MCMC","MTC",priorfile),
                                        usemobility = TRUE,
                                        mobilitydir = "../MobilityMetrics",
                                        isprior = TRUE
)
set.seed(exp(1))
prior_dat <- mcsim(exe_file = exe_file,
                        in_file = prior_template,
                   out_file = file.path(resultsdir,gsub(".in.R",".out",prior_template)),
                 resultsdir = resultsdir)

## Execute: ./../model/mcsim.SEIR.reopen.model.R.exe TX.pred/SEIR_TX_MTC.in.R TX.pred/SEIR_TX_MTC.out

infile_template <- make_infile_template(dat.df,
                                        fips_table,
                                        state_abbr=statenow,
                                        usestatename = TRUE,
                                        createdir = TRUE,
                                        pathdir = resultsdir,
                                        X_iter = "2000",
                                        X_print = "10",
                                        priordir = "../priors",
                                        priorfile = priorfile,
                                        datadatemax = "2020-06-20",
                                        usemobility = TRUE,
                                        mobilitydir = "../MobilityMetrics")
set.seed(exp(1))
out_dat <- mcsim(exe_file = exe_file,
                        in_file = infile_template,
                 resultsdir = resultsdir)#,ignore.stdout = TRUE)

## Execute MCMC: ./../model/mcsim.SEIR.reopen.model.R.exe TX.pred/SEIR_TX_MCMC.in.R TX.pred/SEIR_TX_MCMC1.out

## * Created 'TX.pred/SEIR_TX_MCMC1.out'

## Execute MCMC Check: ./../model/mcsim.SEIR.reopen.model.R.exe TX.pred/SEIR_TX_MCMC.in.R TX.pred/SEIR_TX_MCMC1.out

## * Created 'TX.pred/SEIR_TX_MCMC1.check.out' from the last iteration.

make_diagnostic(out_dat, subset(dat.df,state==statenow), 
                subset(cumdat.df,state==statenow), burnin=0.1,
                pdfname=file.path(resultsdir,
                                  paste0("Test.Prediction.",statenow,".pdf")))

## Warning: Transformation introduced infinite values in continuous x-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 6 rows containing missing values (geom_col).

## Warning: Removed 15 rows containing missing values (geom_path).

## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 1935 rows containing missing values (geom_point).

## quartz_off_screen 
##                 2

statenow <- "TX"
resultsdir <- "TX"
fips_table <- get_fips(datadir = datadir)
dat.df <- get_testdata()
cumdat.df <- get_cumul_testdata()
write.csv(fips_table,
          file=file.path("FIPS_TABLE.csv"),
          row.names = FALSE)
write.csv(dat.df,
          file.path("DAILYTESTDATA.csv"),
          row.names = FALSE)
write.csv(cumdat.df,
          file.path("CUMULTESTDATA.csv"),row.names = FALSE)
burnin <- 0.1
popnow <- fips_table$pop[fips_table$Alpha.code==statenow]
priorfile<-"SEIR.reopen_state_priors_MCMC.in.R"
prior_template <- make_infile_template(dat.df,
                                        fips_table,
                                        state_abbr=statenow,
                                        usestatename = TRUE,
                                        createdir = TRUE,
                                        pathdir = resultsdir,
                                        priordir = "../priors",
                                        priorfile = gsub("MCMC","MTC",priorfile),
                                        usemobility = TRUE,
                                        mobilitydir = "../MobilityMetrics",
                                        isprior = TRUE
)
set.seed(exp(1))
prior_dat <- mcsim(exe_file = exe_file,
                        in_file = prior_template,
                   out_file = file.path(resultsdir,gsub(".in.R",".out",prior_template)),
                 resultsdir = resultsdir)

## Execute: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MTC.in.R TX/SEIR_TX_MTC.out

infile_template <- make_infile_template(dat.df,
                                        fips_table,
                                        state_abbr=statenow,
                                        usestatename = TRUE,
                                        createdir = TRUE,
                                        pathdir = resultsdir,
                                        X_iter = "2000",
                                        X_print = "10",
                                        priordir = "../priors",
                                        priorfile = priorfile,
                                        datadatemax = "2020-06-20",
                                        usemobility = TRUE,
                                        mobilitydir = "../MobilityMetrics")
make_infiles(infile_template,
               exe_file=exe_file,
               chains=1:4,
               useposterior=useposterior,
               resultsdir=resultsdir,
               randomseed=exp(1))

for (chainnum in 1:4) {
  out_dat <- mcsim(exe_file = exe_file,
                          in_file = infile_template,
                   chainnum = chainnum,
                   resultsdir = resultsdir)#,ignore.stdout = TRUE)
}

## Execute MCMC: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC1.out

## * Created 'TX/SEIR_TX_MCMC1.out'

## Execute MCMC Check: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC1.out

## * Created 'TX/SEIR_TX_MCMC1.check.out' from the last iteration.

## Execute MCMC: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC2.out

## * Created 'TX/SEIR_TX_MCMC2.out'

## Execute MCMC Check: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC2.out

## * Created 'TX/SEIR_TX_MCMC2.check.out' from the last iteration.

## Execute MCMC: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC3.out

## * Created 'TX/SEIR_TX_MCMC3.out'

## Execute MCMC Check: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC3.out

## * Created 'TX/SEIR_TX_MCMC3.check.out' from the last iteration.

## Execute MCMC: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC4.out

## * Created 'TX/SEIR_TX_MCMC4.out'

## Execute MCMC Check: ./../model/mcsim.SEIR.reopen.model.R.exe TX/SEIR_TX_MCMC.in.R TX/SEIR_TX_MCMC4.out

## * Created 'TX/SEIR_TX_MCMC4.check.out' from the last iteration.

system(paste("cp",file.path(functiondir,
                            "plot_parameter_results.R"),
             resultsdir))
system(paste("cp",file.path(functiondir,
                            "run_batch_rhat_multicheck.R"),
             resultsdir))
system(paste("cp",file.path(modeldir,basename(exe_file)),
             file.path(resultsdir,
                       gsub(".exe","",basename(exe_file)))))
wd <- getwd()
setwd(resultsdir)
source("plot_parameter_results.R")

##                                [,1]
## GM_NInit.1.                1.278387
## GM_TIsolation.1.           1.007916
## GM_R0.1.                   4.605071
## GM_c0.1.                   1.002404
## GM_TLatent.1.              1.071191
## GM_TRecover.1.             1.757250
## GM_IFR.1.                  1.266413
## GM_TStartTesting.1.        8.136230
## GM_TauTesting.1.           5.839326
## GM_TTestingRate.1.         1.144606
## GM_TContactsTestingRate.1. 1.010410
## GM_TestingCoverage.1.      1.073967
## GM_TestSensitivity.1.      1.008816
## GM_ThetaMin.1.             7.068529
## GM_TauTheta.1.             1.136048
## GM_PwrTheta.1.             1.010689
## GM_HygienePwr.1.           2.895599
## GM_FracTraced.1.           1.230270
## GM_TPosTest.1.             1.123036
## GM_TFatalDeath.1.          5.240737
## GM_TauS.1.                 1.107615
## GM_rMax.1.                 1.025155
## GM_TauR.1.                 1.025190
## alpha_Pos.1.               1.130488
## alpha_Death.1.             1.035886
## LnPrior                    1.852612
## LnData                     4.100665
## LnPosterior                3.277557

## Warning in melt(priors[, 1:(ncol(priors) - 1)], id.vars = 1): The melt generic
## in data.table has been passed a data.frame and will attempt to redirect to the
## relevant reshape2 method; please note that reshape2 is deprecated, and this
## redirection is now deprecated as well. To continue using melt methods from
## reshape2 while both libraries are attached, e.g. melt.list, you can prepend
## the namespace like reshape2::melt(priors[, 1:(ncol(priors) - 1)]). In the next
## version, this warning will become an error.

source("run_batch_rhat_multicheck.R")

##                                [,1]
## GM_NInit.1.                1.278387
## GM_TIsolation.1.           1.007916
## GM_R0.1.                   4.605071
## GM_c0.1.                   1.002404
## GM_TLatent.1.              1.071191
## GM_TRecover.1.             1.757250
## GM_IFR.1.                  1.266413
## GM_TStartTesting.1.        8.136230
## GM_TauTesting.1.           5.839326
## GM_TTestingRate.1.         1.144606
## GM_TContactsTestingRate.1. 1.010410
## GM_TestingCoverage.1.      1.073967
## GM_TestSensitivity.1.      1.008816
## GM_ThetaMin.1.             7.068529
## GM_TauTheta.1.             1.136048
## GM_PwrTheta.1.             1.010689
## GM_HygienePwr.1.           2.895599
## GM_FracTraced.1.           1.230270
## GM_TPosTest.1.             1.123036
## GM_TFatalDeath.1.          5.240737
## GM_TauS.1.                 1.107615
## GM_rMax.1.                 1.025155
## GM_TauR.1.                 1.025190
## alpha_Pos.1.               1.130488
## alpha_Death.1.             1.035886
## LnPrior                    1.852612
## LnData                     4.100665
## LnPosterior                3.277557

setwd(wd)

datadatemax <- "2020-06-20"
folder<-"."
fips_table <- read.csv("FIPS_TABLE.csv",colClasses=c(
  rep("character",4),rep("numeric",2)
))
statenow <- "TX"
scen_model_file<- "SEIR.scenarios.model.R"
scen_exe_file<-makemcsim(scen_model_file,modeldir=modeldir,mdir="../MCSim")

## * Created executable program '../model/mcsim.SEIR.scenarios.model.R.exe'.

output <- run_setpoints1(fips_table,
                           state_abbr=statenow,
                           TPrint=datadatemax,
                           pathdir=folder,
                           scenariosdir = "../scenarios",
                           scenariostemplate=
                             "SEIR.reopen_state_setpoints1_MCMC.in.R",
                           scenarioname = "OneTime",
                           nruns = 0,
                           keepoutfile = TRUE,
                         exe_file=scen_exe_file)

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_OneTime.in ./TX/SEIR_TX_OneTime.out

scen.df <- data.frame(state=sort(rep(statenow,12)),
                      mu_C = rep(rep(c(1,1,2,2),3),length(statenow)),
                      mu_Lambda = rep(rep(c(1,2,1,2),3),length(statenow)),
                      DeltaDelta = rep(c(rep(0,4),rep(0.25,4),rep(-0.25,4)),length(statenow)),
                      stringsAsFactors = FALSE
)
scen.df$scenarioname <- paste("TimeSeries",scen.df$mu_C,scen.df$mu_Lambda,
                              scen.df$DeltaDelta,sep=".")
scen.df$scenariodesc <- paste0(scen.df$mu_C,"X Contact Tracing, ",
                               scen.df$mu_Lambda,"X Testing, ",
                               ifelse(sign(scen.df$DeltaDelta)==1,
                                      paste0("+",100*scen.df$DeltaDelta,"%"),
                                      ifelse(sign(scen.df$DeltaDelta)== -1,
                                      paste0(100*scen.df$DeltaDelta,"%"),
                                      ifelse(sign(scen.df$DeltaDelta)==0,
                                      "Current",""
                                      )))," Reopening")
pdf(file=file.path(folder,statenow,
                   "Scenarios_TestRuns.pdf"),height=4,width=6)
for (j in 1:nrow(scen.df)) {
  scenrow<-scen.df[j,]
  output <- run_setpoints(fips_table,
                           state_abbr=scenrow$state,
                           pathdir=folder,
                          scenariosdir="../scenarios",
                           scenariostemplate=
                             "SEIR.reopen_state_setpoints_MCMC.in.R",
                           scenarioname = scenrow$scenarioname,
                           nruns = 0,
                          mu_C = scenrow$mu_C,
                          mu_Lambda = scenrow$mu_Lambda,
                           DeltaDelta = scenrow$DeltaDelta,
                          rampuptime=14,
                           keepoutfile = FALSE,
                         exe_file=scen_exe_file)
  plot_scenario(alldat.df, output$out_quant,scenrow$state,
                logy=FALSE,
                scenarioname = scenrow$scenariodesc)
}

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.1.1.0.in ./TX/SEIR_TX_TimeSeries.1.1.0.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.1.2.0.in ./TX/SEIR_TX_TimeSeries.1.2.0.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.2.1.0.in ./TX/SEIR_TX_TimeSeries.2.1.0.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.2.2.0.in ./TX/SEIR_TX_TimeSeries.2.2.0.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.1.1.0.25.in ./TX/SEIR_TX_TimeSeries.1.1.0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.1.2.0.25.in ./TX/SEIR_TX_TimeSeries.1.2.0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.2.1.0.25.in ./TX/SEIR_TX_TimeSeries.2.1.0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.2.2.0.25.in ./TX/SEIR_TX_TimeSeries.2.2.0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.1.1.-0.25.in ./TX/SEIR_TX_TimeSeries.1.1.-0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.1.2.-0.25.in ./TX/SEIR_TX_TimeSeries.1.2.-0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.2.1.-0.25.in ./TX/SEIR_TX_TimeSeries.2.1.-0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

## Execute Setpoints: ./../model/mcsim.SEIR.scenarios.model.R.exe ././TX/SEIR_TX_TimeSeries.2.2.-0.25.in ./TX/SEIR_TX_TimeSeries.2.2.-0.25.out

## Warning in melt.data.table(as.data.table(out_dat), id.vars = 1, variable.factor
## = FALSE): 'measure.vars' [GM_NInit, GM_TIsolation, GM_R0, GM_c0, ...] are not
## all of the same type. By order of hierarchy, the molten data value column will
## be of type 'double'. All measure variables not of type 'double' will be coerced
## too. Check DETAILS in ?melt.data.table for more on coercion.

## Warning: Removed 17 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_col).

dev.off()

## quartz_off_screen 
##                 2

## ─ Session info ───────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.1 (2019-07-05)
##  os       macOS Mojave 10.14.6        
##  system   x86_64, darwin15.6.0        
##  ui       X11                         
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Chicago             
##  date     2020-07-03                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────
##  package      * version  date       lib source        
##  assertthat     0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
##  backports      1.1.5    2019-10-02 [1] CRAN (R 3.6.0)
##  bayesplot    * 1.7.1    2019-12-01 [1] CRAN (R 3.6.0)
##  bitops         1.0-6    2013-08-17 [1] CRAN (R 3.6.0)
##  broom          0.5.3    2019-12-14 [1] CRAN (R 3.6.0)
##  callr          3.4.0    2019-12-09 [1] CRAN (R 3.6.0)
##  cellranger     1.1.0    2016-07-27 [1] CRAN (R 3.6.0)
##  cli            2.0.0    2019-12-09 [1] CRAN (R 3.6.0)
##  coda         * 0.19-3   2019-07-05 [1] CRAN (R 3.6.0)
##  colorspace     1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
##  covid19us    * 0.1.3    2020-04-29 [1] CRAN (R 3.6.1)
##  crayon         1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
##  curl           4.3      2019-12-02 [1] CRAN (R 3.6.0)
##  data.table   * 1.12.8   2019-12-09 [1] CRAN (R 3.6.0)
##  DBI            1.1.0    2019-12-15 [1] CRAN (R 3.6.0)
##  dbplyr         1.4.2    2019-06-17 [1] CRAN (R 3.6.0)
##  desc           1.2.0    2018-05-01 [1] CRAN (R 3.6.0)
##  devtools       2.2.1    2019-09-24 [1] CRAN (R 3.6.0)
##  digest         0.6.23   2019-11-23 [1] CRAN (R 3.6.0)
##  dplyr        * 0.8.3    2019-07-04 [1] CRAN (R 3.6.0)
##  ellipsis       0.3.0    2019-09-20 [1] CRAN (R 3.6.0)
##  evaluate       0.14     2019-05-28 [1] CRAN (R 3.6.0)
##  fansi          0.4.0    2018-10-05 [1] CRAN (R 3.6.0)
##  farver         2.0.1    2019-11-13 [1] CRAN (R 3.6.0)
##  forcats      * 0.4.0    2019-02-17 [1] CRAN (R 3.6.0)
##  fs             1.3.1    2019-05-06 [1] CRAN (R 3.6.0)
##  generics       0.0.2    2018-11-29 [1] CRAN (R 3.6.0)
##  GGally       * 1.4.0    2018-05-17 [1] CRAN (R 3.6.0)
##  ggplot2      * 3.2.1    2019-08-10 [1] CRAN (R 3.6.0)
##  ggridges       0.5.1    2018-09-27 [1] CRAN (R 3.6.0)
##  glue           1.3.1    2019-03-12 [1] CRAN (R 3.6.0)
##  gridExtra    * 2.3      2017-09-09 [1] CRAN (R 3.6.0)
##  gtable         0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
##  haven          2.2.0    2019-11-08 [1] CRAN (R 3.6.0)
##  here         * 0.1      2017-05-28 [1] CRAN (R 3.6.0)
##  hms            0.5.2    2019-10-30 [1] CRAN (R 3.6.0)
##  htmltools      0.4.0    2019-10-04 [1] CRAN (R 3.6.0)
##  httr           1.4.1    2019-08-05 [1] CRAN (R 3.6.0)
##  jsonlite     * 1.6      2018-12-07 [1] CRAN (R 3.6.0)
##  knitr          1.26     2019-11-12 [1] CRAN (R 3.6.0)
##  labeling       0.3      2014-08-23 [1] CRAN (R 3.6.0)
##  lattice        0.20-38  2018-11-04 [2] CRAN (R 3.6.1)
##  lazyeval       0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
##  lifecycle      0.1.0    2019-08-01 [1] CRAN (R 3.6.0)
##  lubridate      1.7.4    2018-04-11 [1] CRAN (R 3.6.0)
##  magrittr       1.5      2014-11-22 [1] CRAN (R 3.6.0)
##  memoise        1.1.0    2017-04-21 [1] CRAN (R 3.6.0)
##  modelr         0.1.5    2019-08-08 [1] CRAN (R 3.6.0)
##  munsell        0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
##  nlme           3.1-140  2019-05-12 [2] CRAN (R 3.6.1)
##  pillar         1.4.3    2019-12-20 [1] CRAN (R 3.6.0)
##  pkgbuild       1.0.6    2019-10-09 [1] CRAN (R 3.6.0)
##  pkgconfig      2.0.3    2019-09-22 [1] CRAN (R 3.6.0)
##  pkgload        1.0.2    2018-10-29 [1] CRAN (R 3.6.0)
##  plyr           1.8.5    2019-12-10 [1] CRAN (R 3.6.0)
##  prettyunits    1.0.2    2015-07-13 [1] CRAN (R 3.6.0)
##  processx       3.4.1    2019-07-18 [1] CRAN (R 3.6.0)
##  ps             1.3.0    2018-12-21 [1] CRAN (R 3.6.0)
##  purrr        * 0.3.3    2019-10-18 [1] CRAN (R 3.6.0)
##  R6             2.4.1    2019-11-12 [1] CRAN (R 3.6.0)
##  RColorBrewer   1.1-2    2014-12-07 [1] CRAN (R 3.6.0)
##  Rcpp           1.0.3    2019-11-08 [1] CRAN (R 3.6.0)
##  RCurl        * 1.98-1.2 2020-04-18 [1] CRAN (R 3.6.2)
##  readr        * 1.3.1    2018-12-21 [1] CRAN (R 3.6.0)
##  readxl         1.3.1    2019-03-13 [1] CRAN (R 3.6.0)
##  remotes        2.1.0    2019-06-24 [1] CRAN (R 3.6.0)
##  reprex         0.3.0    2019-05-16 [1] CRAN (R 3.6.0)
##  reshape        0.8.8    2018-10-23 [1] CRAN (R 3.6.0)
##  reshape2     * 1.4.3    2017-12-11 [1] CRAN (R 3.6.0)
##  rlang          0.4.2    2019-11-23 [1] CRAN (R 3.6.0)
##  rmarkdown      2.0      2019-12-12 [1] CRAN (R 3.6.0)
##  rprojroot      1.3-2    2018-01-03 [1] CRAN (R 3.6.0)
##  rstudioapi     0.10     2019-03-19 [1] CRAN (R 3.6.0)
##  rvest        * 0.3.5    2019-11-08 [1] CRAN (R 3.6.0)
##  scales         1.1.0    2019-11-18 [1] CRAN (R 3.6.0)
##  selectr        0.4-2    2019-11-20 [1] CRAN (R 3.6.0)
##  sessioninfo    1.1.1    2018-11-05 [1] CRAN (R 3.6.0)
##  stringi        1.4.3    2019-03-12 [1] CRAN (R 3.6.0)
##  stringr      * 1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
##  testthat       2.3.1    2019-12-01 [1] CRAN (R 3.6.0)
##  tibble       * 2.1.3    2019-06-06 [1] CRAN (R 3.6.0)
##  tidyr        * 1.0.2    2020-01-24 [1] CRAN (R 3.6.0)
##  tidyselect     0.2.5    2018-10-11 [1] CRAN (R 3.6.0)
##  tidyverse    * 1.3.0    2019-11-21 [1] CRAN (R 3.6.0)
##  usethis        1.5.1    2019-07-04 [1] CRAN (R 3.6.0)
##  vctrs          0.2.1    2019-12-17 [1] CRAN (R 3.6.0)
##  viridisLite    0.3.0    2018-02-01 [1] CRAN (R 3.6.0)
##  withr          2.1.2    2018-03-15 [1] CRAN (R 3.6.0)
##  xfun           0.11     2019-11-12 [1] CRAN (R 3.6.0)
##  xml2         * 1.2.2    2019-08-09 [1] CRAN (R 3.6.0)
##  yaml           2.2.0    2018-07-25 [1] CRAN (R 3.6.0)
##  zeallot        0.1.0    2018-01-28 [1] CRAN (R 3.6.0)
## 
## [1] /Users/wchiu/Library/R/3.6/library
## [2] /Library/Frameworks/R.framework/Versions/3.6/Resources/library